/*
 * Copyright (c) 2012, Metron, Inc.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *     * Redistributions of source code must retain the above copyright
 *       notice, this list of conditions and the following disclaimer.
 *     * Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in the
 *       documentation and/or other materials provided with the distribution.
 *     * Neither the name of Metron, Inc. nor the
 *       names of its contributors may be used to endorse or promote products
 *       derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL METRON, INC. BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
package com.metsci.glimpse.util.math.fast;

import com.metsci.glimpse.util.math.LinearInterpolator;

/**
 * Approximates the inverse of <code>erf(x)</code>, valid on the interval
 * [-0.9999779095,0.9999779095], with a maximum error of about <code>6e-5</code>.
 *
 * <p>
 * Being a piecewise linear interpolation of a convex function, the function's
 * magnitude is everywhere an upper bound on the true value, which is fine for
 * the usual application of this function (multipliers for confidence intervals,
 * making the estimate slightly more conservative).
 *
 * @author ellis
 */
public class FastErfInverse
{
    private static final FastErfInverse instance = new FastErfInverse();

    public static final FastErfInverse getInstance()
    {
        return instance;
    }

    private static final double[] XS = { 0, 0.01128341556, 0.02256457469, 0.03384122234, 0.04511110615, 0.05637197780, 0.06762159439, 0.07885771977, 0.09007812584, 0.1012805939, 0.1124629160, 0.1236228962, 0.1347583518, 0.1458671148, 0.1569470331, 0.1679959714, 0.1790118132, 0.1899924612, 0.2009358390, 0.2118398922, 0.2227025892, 0.2335219230, 0.2442959116, 0.2550225996, 0.2657000590, 0.2763263902, 0.2868997232, 0.2974182185, 0.3078800680, 0.3182834959, 0.3286267595, 0.3389081503, 0.3491259948, 0.3592786550, 0.3693645293, 0.3793820536, 0.3893297011, 0.3992059840, 0.4090094534, 0.4187387001, 0.4283923550, 0.4379690902, 0.4474676184, 0.4568866945, 0.4662251153, 0.4754817198, 0.4846553900, 0.4937450509, 0.5027496707, 0.5116682612, 0.5204998778, 0.5292436198, 0.5378986305, 0.5464640969, 0.5549392505, 0.5633233663, 0.5716157638, 0.5798158062, 0.5879229004, 0.5959364972, 0.6038560908, 0.6116812189, 0.6194114619, 0.6270464433, 0.6345858291, 0.6420293274, 0.6493766880, 0.6566277023, 0.6637822027, 0.6708400622, 0.6778011938, 0.6846655502, 0.6914331231, 0.6981039429, 0.7046780779, 0.7111556337, 0.7175367528, 0.7238216140, 0.7300104313, 0.7361034538, 0.7421009647, 0.7480032806, 0.7538107509, 0.7595237569, 0.7651427115, 0.7706680576, 0.7761002683, 0.7814398455, 0.7866873192, 0.7918432468, 0.7969082124, 0.8018828258, 0.8067677215, 0.8115635586, 0.8162710190, 0.8208908073, 0.8254236496, 0.8298702930, 0.8342315043, 0.8385080696, 0.8427007929, 0.8468104962, 0.8508380177, 0.8547842115, 0.8586499465, 0.8624361061, 0.8661435866, 0.8697732972, 0.8733261584, 0.8768031019, 0.8802050696, 0.8835330124, 0.8867878902, 0.8899706704, 0.8930823276, 0.8961238429, 0.8990962029, 0.9020003990, 0.9048374269, 0.9076082860, 0.9103139782, 0.9129555080, 0.9155338810, 0.9180501041, 0.9205051843, 0.9229001283, 0.9252359418, 0.9275136293, 0.9297341930, 0.9318986327, 0.9340079449, 0.9360631228, 0.9380651551, 0.9400150262, 0.9419137153, 0.9437621961, 0.9455614366, 0.9473123980, 0.9490160353, 0.9506732958, 0.9522851198, 0.9538524394, 0.9553761786, 0.9568572531, 0.9582965696, 0.9596950256, 0.9610535095, 0.9623728999, 0.9636540654, 0.9648978648, 0.9661051465, 0.9672767481, 0.9684134969, 0.9695162091, 0.9705856899, 0.9716227333, 0.9726281220, 0.9736026275, 0.9745470093, 0.9754620158, 0.9763483833, 0.9772068366, 0.9780380884, 0.9788428397, 0.9796217795, 0.9803755850, 0.9811049213, 0.9818104416, 0.9824927870, 0.9831525869, 0.9837904586, 0.9844070075, 0.9850028274, 0.9855784998, 0.9861345950, 0.9866716712, 0.9871902752, 0.9876909422, 0.9881741959, 0.9886405487, 0.9890905016, 0.9895245446, 0.9899431565, 0.9903468051, 0.9907359476, 0.9911110301, 0.9914724883, 0.9918207476, 0.9921562228, 0.9924793184, 0.9927904292, 0.9930899398, 0.9933782251, 0.9936556502, 0.9939225709, 0.9941793336, 0.9944262755, 0.9946637246, 0.9948920004, 0.9951114132, 0.9953222650, 0.9955248494, 0.9957194515, 0.9959063484, 0.9960858095, 0.9962580960, 0.9964234618, 0.9965821530, 0.9967344087, 0.9968804605, 0.9970205333, 0.9971548450, 0.9972836068, 0.9974070233, 0.9975252927, 0.9976386070, 0.9977471522, 0.9978511082, 0.9979506491, 0.9980459432, 0.9981371537, 0.9982244380, 0.9983079484, 0.9983878321, 0.9984642313, 0.9985372834, 0.9986071211, 0.9986738725, 0.9987376612, 0.9987986064, 0.9988568234, 0.9989124231, 0.9989655126, 0.9990161950, 0.9990645699, 0.9991107330, 0.9991547766, 0.9991967898, 0.9992368580, 0.9992750637, 0.9993114861, 0.9993462016, 0.9993792835, 0.9994108024, 0.9994408261, 0.9994694199, 0.9994966464, 0.9995225657, 0.9995472358, 0.9995707121, 0.9995930480, 0.9996142945, 0.9996345008, 0.9996537140, 0.9996719792, 0.9996893397, 0.9997058370, 0.9997215109, 0.9997363996, 0.9997505395, 0.9997639656, 0.9997767114, 0.9997888089, 0.9998002889, 0.9998111807, 0.9998215122, 0.9998313105, 0.9998406012, 0.9998494087, 0.9998577566, 0.9998656673, 0.9998731621, 0.9998802615, 0.9998869850, 0.9998933513, 0.9998993781, 0.9999050824, 0.9999104803, 0.9999155873, 0.9999204181, 0.9999249868, 0.9999293067, 0.9999333904, 0.9999372503, 0.9999408977, 0.9999443437, 0.9999475988, 0.9999506730, 0.9999535756, 0.9999563158, 0.9999589021, 0.9999613427, 0.9999636453, 0.9999658172, 0.9999678656, 0.9999697970, 0.9999716177, 0.9999733338, 0.9999749509, 0.9999764744, 0.9999779095 };

    private static final double[] YS = { 0, 0.01000000000, 0.02000000000, 0.03000000000, 0.04000000000, 0.05000000000, 0.06000000000, 0.07000000000, 0.08000000000, 0.09000000000, 0.1000000000, 0.1100000000, 0.1200000000, 0.1300000000, 0.1400000000, 0.1500000000, 0.1600000000, 0.1700000000, 0.1800000000, 0.1900000000, 0.2000000000, 0.2100000000, 0.2200000000, 0.2300000000, 0.2400000000, 0.2500000000, 0.2600000000, 0.2700000000, 0.2800000000, 0.2900000000, 0.3000000000, 0.3100000000, 0.3200000000, 0.3300000000, 0.3400000000, 0.3500000000, 0.3600000000, 0.3700000000, 0.3800000000, 0.3900000000, 0.4000000000, 0.4100000000, 0.4200000000, 0.4300000000, 0.4400000000, 0.4500000000, 0.4600000000, 0.4700000000, 0.4800000000, 0.4900000000, 0.5000000000, 0.5100000000, 0.5200000000, 0.5300000000, 0.5400000000, 0.5500000000, 0.5600000000, 0.5700000000, 0.5800000000, 0.5900000000, 0.6000000000, 0.6100000000, 0.6200000000, 0.6300000000, 0.6400000000, 0.6500000000, 0.6600000000, 0.6700000000, 0.6800000000, 0.6900000000, 0.7000000000, 0.7100000000, 0.7200000000, 0.7300000000, 0.7400000000, 0.7500000000, 0.7600000000, 0.7700000000, 0.7800000000, 0.7900000000, 0.8000000000, 0.8100000000, 0.8200000000, 0.8300000000, 0.8400000000, 0.8500000000, 0.8600000000, 0.8700000000, 0.8800000000, 0.8900000000, 0.9000000000, 0.9100000000, 0.9200000000, 0.9300000000, 0.9400000000, 0.9500000000, 0.9600000000, 0.9700000000, 0.9800000000, 0.9900000000, 1.000000000, 1.010000000, 1.020000000, 1.030000000, 1.040000000, 1.050000000, 1.060000000, 1.070000000, 1.080000000, 1.090000000, 1.100000000, 1.110000000, 1.120000000, 1.130000000, 1.140000000, 1.150000000, 1.160000000, 1.170000000, 1.180000000, 1.190000000, 1.200000000, 1.210000000, 1.220000000, 1.230000000, 1.240000000, 1.250000000, 1.260000000, 1.270000000, 1.280000000, 1.290000000, 1.300000000, 1.310000000, 1.320000000, 1.330000000, 1.340000000, 1.350000000, 1.360000000, 1.370000000, 1.380000000, 1.390000000, 1.400000000, 1.410000000, 1.420000000, 1.430000000, 1.440000000, 1.450000000, 1.460000000, 1.470000000, 1.480000000, 1.490000000, 1.500000000, 1.510000000, 1.520000000, 1.530000000, 1.540000000, 1.550000000, 1.560000000, 1.570000000, 1.580000000, 1.590000000, 1.600000000, 1.610000000, 1.620000000, 1.630000000, 1.640000000, 1.650000000, 1.660000000, 1.670000000, 1.680000000, 1.690000000, 1.700000000, 1.710000000, 1.720000000, 1.730000000, 1.740000000, 1.750000000, 1.760000000, 1.770000000, 1.780000000, 1.790000000, 1.800000000, 1.810000000, 1.820000000, 1.830000000, 1.840000000, 1.850000000, 1.860000000, 1.870000000, 1.880000000, 1.890000000, 1.900000000, 1.910000000, 1.920000000, 1.930000000, 1.940000000, 1.950000000, 1.960000000, 1.970000000, 1.980000000, 1.990000000, 2.000000000, 2.010000000, 2.020000000, 2.030000000, 2.040000000, 2.050000000, 2.060000000, 2.070000000, 2.080000000, 2.090000000, 2.100000000, 2.110000000, 2.120000000, 2.130000000, 2.140000000, 2.150000000, 2.160000000, 2.170000000, 2.180000000, 2.190000000, 2.200000000, 2.210000000, 2.220000000, 2.230000000, 2.240000000, 2.250000000, 2.260000000, 2.270000000, 2.280000000, 2.290000000, 2.300000000, 2.310000000, 2.320000000, 2.330000000, 2.340000000, 2.350000000, 2.360000000, 2.370000000, 2.380000000, 2.390000000, 2.400000000, 2.410000000, 2.420000000, 2.430000000, 2.440000000, 2.450000000, 2.460000000, 2.470000000, 2.480000000, 2.490000000, 2.500000000, 2.510000000, 2.520000000, 2.530000000, 2.540000000, 2.550000000, 2.560000000, 2.570000000, 2.580000000, 2.590000000, 2.600000000, 2.610000000, 2.620000000, 2.630000000, 2.640000000, 2.650000000, 2.660000000, 2.670000000, 2.680000000, 2.690000000, 2.700000000, 2.710000000, 2.720000000, 2.730000000, 2.740000000, 2.750000000, 2.760000000, 2.770000000, 2.780000000, 2.790000000, 2.800000000, 2.810000000, 2.820000000, 2.830000000, 2.840000000, 2.850000000, 2.860000000, 2.870000000, 2.880000000, 2.890000000, 2.900000000, 2.910000000, 2.920000000, 2.930000000, 2.940000000, 2.950000000, 2.960000000, 2.970000000, 2.980000000, 2.990000000, 3.000000000 };

    private static final LinearInterpolator INTERPOLATOR = LinearInterpolator.createFromArrays( XS, YS );

    private FastErfInverse()
    {
    }

    public double evaluate( double x )
    {
        if ( x >= 0.0 ) return INTERPOLATOR.evaluate( x );
        return -INTERPOLATOR.evaluate( -x );
    }
}
